rm(list=ls(all=TRUE)) # clear workspace

require(rms) # for logit models
require(ggplot2) # for plotting
require(survival) # for conditional logit

setwd("D:/Dropbox/Reeder and Smith/Revisions/Final Documents/Replication") # set working directory
ann <- read.csv("ReederSmithFPA_annualcounts.csv") # bring in annual counts for Figure 1
core <- read.csv("ReederSmithFPA_rep.csv") # bring in data for analysis

# prep dist for core
ddist <- datadist(core)
options(datadist='ddist')

# Create Figure 1 
figure1a <- ggplot() + geom_area(aes(y = strikes, x = year), data = ann, stat="identity", fill="#999999") +
  theme_minimal() + labs(x="Year", y="US Strikes")
figure1a 
figure1b <- ggplot() + geom_area(aes(y = killings, x = year), data = ann, stat="identity", fill="#999999") +
  theme_minimal() + labs(x="Year", y="Targeted Civilian Killings by Al-Shabaab")
figure1b 

# Table 1 models
lr.agg <- lrm(tarkil_dum~strike_lag+sommil_lag+amisom_lag+kenya_lag+ethiopia_lag+capture_lag+ter_lag+pcewks+ttime_mean+borders+capital, data=core, x=T, y=T)
lr.agg <- robcov(lr.agg, cluster=core$gridid) # cluster errors
lr.agg # model 1 results
clogit(tarkil_dum~strike_lag+sommil_lag+amisom_lag+kenya_lag+ethiopia_lag+capture_lag+ter_lag+pcewks+strata(gridid), data=core) # model 2 results
summary(lm(tarkil_dum~strike_lag+sommil_lag+amisom_lag+kenya_lag+ethiopia_lag+capture_lag+ter_lag+pcewks+factor(gridid), data=core)) # model 3 results
summary(lm(tarkil_dum~strike_lag+sommil_lag+amisom_lag+kenya_lag+ethiopia_lag+capture_lag+ter_lag+pcewks+factor(gridid)+factor(yr)+factor(mo)+factor(wk), data=core)) # model 4 results

# Capture odds ratios and create Figure 3 (using model 1)
eff.df <- as.data.frame(summary(lr.agg))
eff.df <- eff.df[eff.df$Type == 2,]
eff.df$variable <- c("US Strikes (t-1)", "Somalia Battles (t-1)", "AMISOM Battles (t-1)", "Kenya Battles (t-1)", "Ethiopia Battles (t-1)",
                     "Resource Capture (t-1)", "Territorial Loss (t-1)", "Peace Weeks", "Urban Center", "Distance to Border", "Distance to Capital")
eff.df <- eff.df[c("variable", "Effect", "Lower 0.95", "Upper 0.95")]
names(eff.df) <- c("variable", "or", "low", "high")
eff.df$value <- round(eff.df$or, digits=3)
eff.df$sig <- "**"
eff.df$sig[eff.df$low < 1 & eff.df$high > 1] <- ""
eff.df$label <- paste(eff.df$value, eff.df$sig, sep="")
eff.df$variable <- factor(eff.df$variable, levels=rev(eff.df$variable))
#
figure3 <- ggplot(data=eff.df, aes(x=variable, y=or, ymin=low, ymax=high)) +
           geom_pointrange() + geom_text(aes(label=label), vjust="bottom", nudge_x = 0.1) +
           geom_hline(yintercept=1, lty=2, colour="black") +  # add a dotted line at x=1 after flip
           coord_flip() +  # flip coordinates (puts labels on y axis)
           xlab("") + ylab("Odds Ratios (95% CI)") +
           theme_minimal()  # use a white background
figure3

# Table 2 models
lr.dis <- lrm(tarkil_dum~assets_lag+decap_lag+memb_lag+sommil_lag+amisom_lag+kenya_lag+ethiopia_lag+capture_lag+ter_lag+pcewks+ttime_mean+borders+capital, data=core, x=T, y=T)
lr.dis <- robcov(lr.dis, cluster=core$gridid) # cluster errors
lr.dis # model 5 results
clogit(tarkil_dum~assets_lag+decap_lag+memb_lag+sommil_lag+amisom_lag+kenya_lag+ethiopia_lag+capture_lag+ter_lag+pcewks+strata(gridid), data=core) # model 6 results
summary(lm(tarkil_dum~assets_lag+decap_lag+memb_lag+sommil_lag+amisom_lag+kenya_lag+ethiopia_lag+capture_lag+ter_lag+pcewks+factor(gridid), data=core)) # model 7 results
summary(lm(tarkil_dum~assets_lag+decap_lag+memb_lag+sommil_lag+amisom_lag+kenya_lag+ethiopia_lag+capture_lag+ter_lag+pcewks+factor(gridid)+factor(yr)+factor(mo)+factor(wk), data=core)) # model 8 results

# Capture odds ratios and create Figure 4 (using model 5)
eff2.df <- as.data.frame(summary(lr.dis))
eff2.df <- eff2.df[eff2.df$Type == 2,]
eff2.df$variable <- c("Asset Strikes (t-1)", "Decapitation Strike (t-1)", "Membership Strike (t-1)", "Somalia Battles (t-1)", "AMISOM Battles (t-1)", "Kenya Battles (t-1)", "Ethiopia Battles (t-1)",
                      "Resource Capture (t-1)", "Territorial Loss (t-1)", "Peace Weeks", "Urban Center", "Distance to Border", "Distance to Capital")
eff2.df <- eff2.df[c("variable", "Effect", "Lower 0.95", "Upper 0.95")]
names(eff2.df) <- c("variable", "or", "low", "high")
eff2.df <- eff2.df[1:3,]
eff2.df$value <- round(eff2.df$or, digits=3)
eff2.df$sig <- "**"
eff2.df$sig[eff2.df$low < 1 & eff2.df$high > 1] <- ""
eff2.df$label <- paste(eff2.df$value, eff2.df$sig, sep="")
eff2.df$variable <- factor(eff2.df$variable, levels=rev(eff2.df$variable))
#
figure4 <- ggplot(data=eff2.df, aes(x=variable, y=or, ymin=low, ymax=high)) +
           geom_pointrange() + geom_text(aes(label=label), vjust="bottom", nudge_x = 0.05) +
           geom_hline(yintercept=1, lty=2) +  # add a dotted line at x=1 after flip
           coord_flip() +  # flip coordinates (puts labels on y axis)
           xlab("") + ylab("Odds Ratios (95% CI)") +
           theme_minimal()  # use a white background
figure4


